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AN IMPLICIT ALGORITHM FOR THE CONSERVATIVE, TRANSONIC FULL POTENTIAL 


EQUATION WITH EFFECTIVE ROTATED DIFFERENCING 
Terry L. Holst 
Ames Research Center 
and 

John Albert* 

University of Santa Clara 
SUMMARY 


A new differencing scheme for the conservative full potential equation 
which effectively simulates rotated differencing is presented. The scheme is 
implemented by an appropriate upwind bias of the density coefficient along 
both coordinate directions. A fast, fully implicit, approximate factorization 
iteration scheme is then used to solve the resulting dxfference equations. 
Solutions for a number of traditionally difficult transonic airfoil test cases 
are presented. 


I. INTRODUCTION 


In reference 1, a new implicit approximate factorization algorithm (AF2) 
was presented. This algorithm has been applied to the solution of the tran- 
sonic small-disturbance equation (ref. 2) and the conservative full potential 
equation (refs. 3 and 4). In both cases, significant improvements in conver- 
gence speed have been realized over standard successive line overrelaxation 
'SLOR) algorithms. Stability in the full potential formulations for super- 
sonic regions of flow has been achieved by the addition of an artificial vis- 
cosity term similar to that introduced in reference 5. However, in the pres- 
ent formulation the addition of the artificial viscosity term has been achieved 
by an upwind bias of the density coefficient. This strategy greatly simplifies 
the solution procedure and effectively allows the simple two- and three-banded 
matrix form of the AF scheme to be retained over the entire flow field, even 
in regions of supersonic flow. Other researchers (refs. 5-8) have used simi- 
lar steady-state differencing procedures in a wile variety of problems to fur- 
ther substantiate this approach as being both “ciiiable and flexible. 

The use of a numerical transformation to establish an arbitrary finite- 
difference mesh Woo introduced as an aspect of the present full potential 
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algorithm in reference 4. This grid generation technique was developed in 
reference 9 and has been extensively used for many applications (refs. 10 
and 11). The basic transformation (x,y -> ^,n) is illustrated in figure 1. 

With this approach, body-fitted finite-difference meshes for arbitrary airfoil 
geometries can be routinely and efficiently generated. Extension of this pro- 
cedure to more complicated geometries including three-dimensional problems is 
also possible. 

In reference 4, the density coefficient was upwinded along only the wrap- 
around direction coordinate, see fig. 1) and not in the n-coordinate 
direction. Upwinding the density coefficient along only one coordinate direc- 
tion does not in every case allow the numerical domain of dependence to fully 
accommodate the physical domain of dependence. Hence, in certain cases (usu- 
ally solutions with strong shock waves), convergence is difficult or impossible 
to achieve. The purpose of this investigation is to develop an effective form 
of "rotated" differencing which will allow the numerical domain of dependence 
to more nearly match the physical domain of dependence, and thereby, improve 
the reliability of the present algorithm. 

The finite-difference scheme including the new rotated differencing is 
presented in the next section along with several helpful guidelines for improv- 
ing the convergence speed of the AF2 algorithm. The last section contains a 
variety of computational results including several classically difficult test 
cases which subject the new differencing scheme to the extreme test. 

The authors express their gratitude to Vladimir Drobot, University of 
Santa Clara, for his assistance during the course of this study. 


II. THE FULL POTENTIAL EQUATION ALGORITHM 


A. Governing Equations 

The full potential equation written in strong conservation-law form is 
given by 


(p({) ) + (p 4> ) =0 (la) 

p = [i - (♦/ + ♦/)] (lb) 

where the density (p) and velocity components (({ij^ and cj)y) are nondimensional- 
ized by the stagnation density (Pg) and the critical sound speed (a^) , respec- 
tively; X and y are Cartesian coordinates; and y is the ratio of specific 

heats . 

Equation (1) expresses mass conservation for flows that are isentropic 
and irrotatioual. The corresponding shock-jump conditions are valid 
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approximations to the Ranklne-Hugonlot relations for many transonic flow 
applications. A comparison of the Isentroplc and Ranklne-Hugonlot shock polars 
Is given In reference 12 . 


Equation (1) Is transformed from the physical domain (Cartesian coordi- 
nates) Into a computational domain by using a general Independent variable 
transformation. This general transformation. Indicated by (see fig. 1) 

^ = Ux,y) \ 

I 

n = n(x,y) ) 

maintains the strong conservation-law form of equation ( 1 ) as discussed In 
references 10 and 13-15. The full potential equation written In the computa- 
tional domain (^-n coordinate system) Is given by 



(3a) 


where 


U = + A 2 (j)^ 

V = A^(j)^ + A 3 <|)^ 

Ai = C ^ + C ^ 

^ XX y y 

A3 = n 2 + ^ 2 

^ X y 


(3b) 


(4) 


U and V are the contravarlant velocity components along the ^ and n direc- 
tions, respectively; A^ , A 2 , and A 3 are metric quantities; and J Is the 
Jacobian of the transformation. 


The transformed full potential equation (eq. (3)) Is only slightly more 
complicated than the original Cartesian form (eq. (1)) and offers several sig- 
nificant advantages. The main advantage Is that boundaries associated with 
the physical domain are transformed to boundaries of the computational domain. 
This aspect Is Illustrated In figure 1, where the physical and computational 
domains for a typical transformation are shown. The Inner airfoil boundary 
becomes the q = nmax computational boundary. Note that no restrictions have 
been placed on the shape of the outer boundary. Arbitrarily shaped outer 
boundaries. Including wind-tunnel walls, may be used. 
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B. Grid Generation 


The grid generation algorithm used in the present study is an adaptation 
of the scheme introduced by Thompson, Thames, and Mastin (ref. 9) and has been 
discussed in detail in reference A. Basically, this scheme uses numerically 
generated solutions of Laplace's equation (or in some cases Poisson's equation) 
to establish regular and smooth finite-difference meshes around arbitrary 
bodies. These equations are transformed to (and solved in) the computational 
domain (i.e., ^ and n are the independent variables and x and y are the 
dependent variables). A fast approximate factorization relaxation algorithm 
is used to solve these equations. Once values of x and y are known for each 
point in the mesh, the metric quantities of equation (A) are computed by stan- 
dard finite-difference formulas. Further details about this procedure are 
found in reference A. 


C. Spatial Differencing 


A second-order accurate finite-difference approximation to the full 
potential equation (eq. (3a)) is given by 





+ 



l.j+l/2 


0 


(5a) 


where 


i+l/2,j 


^ ^^i+l/2,j . j+1 ‘^i+l,j-l 




(5b) 




‘*’i+i,j " “^i-i,: 


(5c) 


The quantities p, Aj , A 2 , A 3 , and J are all stored at integer points in the 
finite-difference mesh (i.e., at i,j). Values needed at half points (i.e., 
i+l/2,j or i,j+l/2) are obtained by using simple averages. The operators, 

6 - ( ) and 6 ( ) 

A n 


A 



are first-order accurate backward difference operators in the ^ and n direc- 
tions, respectively, and are defined by 


>ij ■ < >i,j - ‘ 


n i,j i,j i.j-1 


The density calculation is performed in a straightforward manner by using equa- 
tion (3b). Values of U, V, i})^, and (fi^ required for the density calculation 
are given by 



^i,j ■ ^ 

Vi,j = (A2<^^ + 


( 6 ) 


Special formulas replace equations (5) and (6) at boundaries and are discussed 
in reference 4. 


Equation (5) is a suitable finite-difference scheme for subsonic flow 
regions. However, for supersonic regions, a properly chosen artificial vis- 
cosity term must be added. For example, Jameson (ref. 5) adds the following 
viscosity term to the Cartesian form of the full potential equation: 

-Ax(vp (|) ) - Ay(vp |(j) I) (7) 

X X X y y y 

where v = max [0,1 - (1/M^)] (M is the local Mach number). Exact implementa- 
tion of equation (7) in the present case involving the ^ - n coordinate 
system is difficult. Approximate implementation was achieved in reference 4 
by using an artificial viscosity of the following form 

-i5(vpjl2l)^ (8) 

where the absolute value of the U-velocity component has been included as a 
consequence of the wraparound coordinate system. Equation (8) is a good 
approximation to the first term of equation (7) on the upper and lower sur- 
faces of an airfoil where the general 5 - n coordinate system is approxi- 
mately Cartesian. This strategy has been very successful for cases in which 
the supersonic zones are relatively small (i.e., cases in which the leading 
and trailing edges are reasonably far removed from supersonic flow). However, 
for cases containing strong shock waves at the trailing edge or supersonic 
free streams, a more complete approximation to equation (7) is required. An 
obvious extension of equation (8) is given by 
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( 9 ) 


As pointed out by Jameson (ref. 5), addition of equation (7) (with appropriate 
upwind differencing) to the centrally differenced full potential equation 
(Cartesian form) provides automatic upwind differencing of the streamwise 
terms in supersonic regions, thus including the full effect of rotated differ- 
encing. Use of equation (9) (with appropriate upwind differencing), there- 
fore, closely approximates the effects of a rotated differencing scheme. 

The complete finite-difference approximation to equation (3a) is given by 


(~) 


i+l/2,j 


■ ■ ‘'i,j+2+l/2j ■ ° 


( 10 ) 


where 


k = 


i = 


’ -1 

when 

^i+l/2,j 

, 1 

when 


-1 

when 

'^i,j+l/2 

, 1 

when 

j+ 1/2 


> 0 

< 0 
> 0 

< 0 


(11a) 


(11b) 


In the artificial viscosity term, is evaluated with a backward difference 

when j ^ ® v’lth a forward difference when ^2 j ^ Like- 
wise, is’evaluated with a backward difference when Vi^j+l/2 ^ 0 and with 

a forward difference when 0. This maintains an upwind influence 

in the differencing scheme for supersonic regions anywhere in the finite- 
difference mesh for any orientation of the velocity vector. The test on 
^i+l/2,j is determined by whether a grid point is in the upper or lower half 
of tne’airfoil flow field, thus simplifying the computational algorithm. 


The scheme given by equations (10) and (11) is centrally differenced and 
second-order accurate in subsonic regions. In supersonic regions, the differ- 
encing is a combination of the second-order accurate central differencing used 
in subsonic regions and the first-order accurate upwind differencing resulting 
from the addition of artificial viscosit/. As the flow becomes increasingly 
supersonic, the scheme is Increasingly retarded in the upwind direction. 
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As shown in reference 3, equation (10) can be rearranged to give 



where 

^i “ ^i+l/2,j^i+k+i/2,j 

Pj = [(1 - + ^i,j+i/2^i,j+il+l/2 


(12a) 


(12b) 

(12c) 


The addition of the artificial viscosity given by equation (9) is thus equiva- 
lent to retarding the density in equation (5a). Artificial viscosity is not 
added explicitly as in the Jameson procedure. 

In the n direction near the airfoil boundary, a problem arises when the 
V-component of velocity is negative. The upwind biased density coefficient at 
i,NJ-l/2 (pj^j_i) requires the use of the density at i,NJ+l/2, which is out- 
side the finite difference mesh (see eq. (12c)). A value of density at this 
location is extrapolated from interior values using 

^i,NJ+l/2 " " ^°^i,NJ-l ^'^i,NJ-2^^® 


Other values of density exterior to the finite-difference mesh are not 
required since at the airfoil surface (j = NJ) the following boundary condi- 
tion is applied: 


^NJ-1 



i,NJ-l/2 


-P 



i,NJ+l/2 


(13a) 


This boundary condition applied in either subsonic or supersonic regions is 
first-order accurate. Another boundary condition, which is second-order accu- 
rate for subsonic flow and first-order accurate for supersonic flow, has been 
tried and can be expressed as 


(“) 


8 _ /V\ 

3 ^NJ-l/2\j/ 


- 3p 


i,NJ 


i,NJ 


NJ-1 \J/. 


i,NJ-l/2 


I ^NJ 


-S) 


l,NJ-3/2 

(13b) 


where the first term is zero because of the = 0 boundary condition. 

Results using these two boundary conditions have been compared and displayed 
little or no difference. 

Several variations of the spatial difference scheme given by equation (12) 
have been considered. These variations arise from different evaluations of 
the artificial viscosity coefficient v and are presented as follows: 
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Version 1 




''i+l/2,j ~ ^”i+l/2,j " 


(14a) 


”i+l/2.j " 


^i+l/2,j 


= 0 


(14b) 


Version 2 


"’i+l/2,j ' 


.•naxKM^ j - 1)C,0] for j ^ ^ 


(15) 


Uax[(M 2 ^^^^ - DC.O] for < 0 


where the standard supersonic form for v (i.e., v = 1 - 1/M^) has been multi- 
plied by CM^. The parameter C is a user-specified constant, usually 
between 1.0 and 2.0. Multiplication by this additional term is used to 
increase the amount of artificial viscosity and therefore, the amount of 
upwinding in the difference scheme. In addition, the value of v is never 
allowed to exceed 1.0. This seems to improve the stability and in some cases 
improves the convergence rate. Expressions for '^i,j+i /2 are written simi- 
larly to those for \'i+i/ 2 ,j and depend on Mi,j+i /2 and 

In addition to the basic form of artificial viscosity already presented 
(see eq. (9)), another form was tested and is given by 


-A 5 (vpj ^ * 5 )^ - Ai(vp^ ^ (16) 

The complete finite-difference approximation to equation (3a) using this form 
of artificial viscosity is obtained similarly to equation ( 10 ) and will not be 
presented. The coefficient v for this artificial viscosity form is evalu- 
ated similarly to the v coefficient of version 2 (see eq. (15)). The only 
difference is that the sign test on V^^ j+l /2 replaced by a sign test on 
(‘t*r))i j+ 1 / 2 ’ is desirable (for this case) because both the sign of 41 ^ 

and the density upwind direction should change at the same point to maintain 
consistency. This artificial viscosity form along with version 2 evaluation 
of V will be referred to as version 3 in the section on computed results. 


D. The AF2 Iteration Scheme 

The AF2 fully implicit approximate factorization scheme used in the pres- 
ent study is discussed in reference 4. Implementation of the scheme is 
achieved by writing it in a two-step form given by 
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i 


Step 1 


Step 2 



( 17 ) 


(18) 


where the n superscript is an iteration index, u) is a relaxation parameter 
(set equal to 1.8 for all cases), is the nth iteration residual oper- 

ator (defined by eq. (12a)), and f^ j is an intermediate result stored at 
each point in the finite-difference mesh. In step 1, the f array is obtained 
by solving a simple bidiagonal matrix equation for each C = constant line. 

The correction array (C^.j = 4>it^i ~ *^i 1^ then obtained in the second step 
from the f array by solving a ’tridia^onal matrix equation for each 
n = constant line. Note that with the AF2 scheme the q-direction differ- 
ence approximation is split between the two steps. This generates a 4>nt iyP® 
term, which is useful to the iteration scheme as timelike dissipation. (The 
iterative process is considered as an iteration in pseudotime. Thus, the time 
derivative is introduced by ( - ( )*' ~ ( )^ . ) The split q term also 

places a sweep direction restriction on both steps, namely, outward (away from 
the airfoil) for the first step and inward (toward the airfoil) for the second 
c,tep. No sweep restrictions are placed on either of the two sweeps due to 
flow direction. 


A type term has been added inside the brackets of step 2 (see 

eq. (18)), to provide time-dependent dissipation in the ^ direction. The 
parameter 3 is determined as follows 


/ 



^M>1 

^M<1 


if 


if 


M. . < 1 

1-1. J 

< 1 

1+1 ,j 


upper surface 
lower surface 
upper surface 
lower surface 


(19) 


where 3j^<^ is fixed at a value of 0.3 and 6j^>i is a user specified con- 
stant which can be adjusted as needed. (The default value for B^m Is equal 
to 1.0, which is sufficient to stabilize most moderate strength shock-wave 
calculations.) The double arrow notation on the ^-difference operator indi- 
cates that the difference is always upwind, which on the upper surface is a 
backward difference and on the lower surface is a forward difference. The 
sign is chosen in such a way that the addition of increases the magni- 

tude of the second sweep diagonal. 


The quantity a appearing in equations (17) and (18) can be considered 
as At“^. This direct analogy to time provides one strategy for obtaining 
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fast convergence, namely, advance time as fast as possible with large time 
steps (i.e., small values of a). As pointed out in reference 3, this is 
effective for attacking the low-frequency errors but not the high-frequency 
errors. The best overall approach is to use an a sequence containing sev- 
eral values of a. The small values are particularly effective for reducing 
the low-frequency errors, and the large values are particularly effective for 
reducing the high-frequency errors. A suitable a sequence with analytically 
estimated endpoints is given in reference 3. 

An additional possibility for Improving the convergence rate is to not 
only use an a sequence, but also vary the a sequence with spatial porition. 
With such a strategy, the a sequence can be adjusted locally to accommodate 
any type of error, i.e., a relatively large residual would indicate a high- 
frequency error and the need for smaller alphas, and a relatively small resid- 
ual would indicate low-frequency error and the need for larger alphas. To 
date, this strategy has not been implemented in any sophisticated manner. 
However, a simplified version has been tested which assumes that the highest 
frequency errors will be supported in the finest region of the finite- 
difference mesh (i.e., near the airfoil). In this region (specifically for 
the three n = constant coordinate lines nearest the airfoil), the lowest 
alphas (largest time steps) are never alio- ed to fall below a certain value. 
With this constraint, it was found that the lowest global values of alpha 
could be reduced by as much as a factor of five which resulted in a convergence 
rate improvement of roughly a factor of two. 

Normally, flow-field, type-dependent differencing is used to achieve sta- 
bility in transonic flow calculations. Incorporating these different operators 
into iteration procedures, such as the AF2 scheme presented here, would be 
cumbersome if not impossible. Using the upwind bias of the density coeffi- 
cient, which is always evaluated at the nth iteration level, allows the simple 
two- and three-banded matrix form of the AF2 scheme to be retained over the 
entire flow field, even in regions of supersonic flow. In fact, use of 
upwinded density coefficients in any general iteration scheme could be used to 
remove the difficulties introduced by type-dependent differencing. The result- 
ing general scheme would retain the same basic differencing (at the n + 1 
iteration level) throughout the entire flow field, relying on the upwind bias 
of the density (at the nth iteration level) to provide the artificial viscosity 
in supersonic flow regions. This represents a significant simplification in 
the handling of supersonic flow regions for transonic flow calculations. 


III. COMPUTED RESULTS 


The new rotated differencing algorithm introduced in the previous section 
is evaluated in this section by presenting a range of numerically computed 
examples, including transonic and supersonic free-stream test cases. In all 
cases the numerically generated finite-difference mesh was converged until the 
maximum residual dropped by three orders of magnitude. For a typical mesh 
with 4470 grid points (149x30), this required approximately 30 to 40 Iterations 
and about 5 sec of CPU time on the Ames 7600 computer. An example of a 
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numerically generated finite-difference mesb for the NACA 0012 airfoil is 
shown in figure 2. 


A. Artificial Viscosity Variations 

A series of calculations has been computed to determine the effect of the 
spatial differencing variations presented in the previous section. The test 
case chosen for this purpose is an NACA 0012 airfoil at a free-stream Mach 
number of 0.75 and 2 ° angle of attack. Systematic variations of the param- 
eter C for all three artificial viscosity versions (see eqs. (14)-(16)), 
including both rotated (NDIF = 1) and nonrotated (NDIF = 0) forms, were inves- 
tigated. (The nonrotated difference scheme simply has no n-direction arti- 
ficial viscosity term. This causes the term (pV/J)^ to remain centrally dif- 
ferenced in supersonic regions.) These solutions are typically represented by 
the pressure coefficient distributions shown in figures 3 and 4. As expected, 
cases involving larger amounts of artificial viscosity (i.e., larger values of 
C) produce slightly smaller lift coefficients and larger amounts of shock 
smearing. Use of the rotated differencing option (NDIF = 1) had a similar 
effect causing a slight amount of additional shock smearing. This trend was 
basically exhibited by each of the three differencing variations. Two of the 
schemes (versions 2 and 3) produced similar solutions for most of the cases. 

The version 1 solutions (see fig. 3) were somewhat different, displaying dis- 
persive overshoots at the shock wave. Increasing the value of C reduces the 
size of the overshoot but does not entirely eliminate it. Some difficulty in 
maintaining stability as well as the normal convergence rate has been experi- 
enced with the version 1 scheme and is attributed to the existence of disper- 
sion in these solutions. 

The lack of an overshoot for the version 2 solutions (see fig. 4) can be 
explained as follows: The artificial viscosity coefficient /2 j ^s com- 

puted in the version 2 differencing scheme is actually evaluated at i,j (for 
Ui+i/ 2 ,j > 0) instead of at i+l/2,j. This effectively provides a substantial 
increase in the artificial viscosity at the shock wave, especially at the last 
supersonic point near the shock-wave center (because j^’^i+l/2 j ^ 
this point), and usually a slight decrease in other supersonic regions 
(because ’^i , j /^i+l /2 , j 1) • The effect of this substantial increase in 

artificial viscosity at the shock is to dramatically inhibit the formation of 
dispersive shock-wave overshoots. Since the version 1 differencing scheme 
does not have this property, shock-wave overshoots can be suppressed only by 
adding a large amount of artificial viscosity everywhere in the supersonic 
region. 

Numerical experiments with other cases including several classically dif- 
ficult cases, which will be oresented in the next section, have been conducted 
with artificial viscosity versions 2 and 3. Although the differences were 
small, the version 2 differencing was chosen as superior. The main advantage 
was slightly faster and smoother convergence especially for supersonic free- 
stream cases. From this point on only artificial viscosity version 2 solutions 
will be presented. 
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B. Convergence Reliability 


Three clasically difficult test cases have been investigated to help aval- 
uate the convergence reliability of the present algorithm, especially the 
newly implemented rotated differencing option. These cases involve three dif- 
ferent free-stream conditions for an NACA 0012 airfoil: (1) = 0.98, at = 0“, 

(2) Moo = 0.95, a = 4°, and (3) Moo “ 1.15, a = 0°. Each case has an unusually 
large amount of supersonic flow and a rather unique shock-wave pattern. Clas- 
sical relaxation schemes would probably experience either very slow convergence 
or instability with these cases. 

Solutions with approaching unity— The pressure coefficient distribu- 

tion for the first case (Moo = 0.98, a = 0°) is shown in figure 5. This case 
is nonlifting (Cj^ = 0.0 and Cj^ = 0.0) and has a large amount of wave drag 
(Cp = 0.1038). The airfoil is almost totally immersed in supersonic flow 
which begins at approximately 5% of chord and continues through an oblique 
shock wave at the trailing edge. A more complete description of the flow 
field is shown by the Mach number contour map displayed in figure 6. Existence 
of a rather unique double shock-wave stiucture, sometimes referred to as a 
"fishtail shock," is clearly evident. Relatively weak oblique shocks emanate 
from the trailing edge and merge with a normal shock downstream of the airfoil. 
The triangular region between the oblique and normal shocks has a nearly con- 
stant supersonic Mach number (approximately equal to 1.1). This shock-wave 
pattern is characteristic of solutions with free-stream Mach numbers near 
unity. It has been observed experimentally as well as computationally and is 
generally considered to be the correct qualitative solution. For instance, a 
fishtail shock solution, for a 10% circular arc airfoil at M^,, = 0.98 and 
zero degrees angle of attack, was presented in reference 6. This calculation 
was a solution to the conservative full potential equation using a Cartesian 
mesh and small-disturbance boundary conditions. Because the flow was essen- 
tially aligned with the finite-difference mesh, "rotated differencing" was not 
necessary. However, wi*.'i the present wraparound coordinate system, rotated 
differencing is essential for maintaining the stability of difficult cases such 
as the present fishtail shock solution. For example, in the supersonic trian- 
gle between the oblique and normal shocks (see fig. 6), the 5 direction is 
nearly normal to the direction of flow. This makes the ^-direction artifi- 
cial viscosity term ineffective in achieving an upwind influence. The upwind 
bias of the density along the n coordinate (which is achieved by a forward 
difference of in the q-direction artificial viscosity term) provides 

the necessary upwind influence to stabilize the calculation. 

Another very important Ingredient for maintaining stable convergence is 
having the proper amount of timelike dissipation (<t'^f and <J>r^) • (This fact 
concerning AF schemes was first mentioned to the authors by J. South, NASA 
Langley.) A fixed amount of is automatically added to the iterative 

process by the construction of the present AF2 algorithm. Likewise, the 
amount of added in subsonic regions is also fixed = 0.3, see 

eq. (19)). The amount of added in supersonic regions is controlled by 

a user-specified constant (6j.j>j) which can be adjusted as needed. For the 
fishtail shock calculation just presented, the parameter was equal to 

five. Smaller values produced instabilities for this calculation. 
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As pointed out in reference 16, nonconservative full potential equation 
solutions with the free-streara Mach number approaching unity are characterized 
by strong oblique shock waves at the trailing edge followed by subsonic flow. 
The fishtail shock structure for these cases is not predicted. Conservative 
vs nonconservative differencing was the subject of discussion in reference 17 
where similar differences for the transonic small-disturbance equation were 
reported. It is generally understood that these differences are the result of 
effective mass creation at shock waves for the nonconservative differencing 
schemes. Therefore, to obtain the proper mass balance and the correspondingly 
correct solution (to either the full potential or transonic small-disturbance 
potential equation), conservative differencing is required. 

The pressure coefficient distribution for the second case (also involving 
a tree-stream Mach number near unity (Moo = 0.95)) is shown in figure 7. This 
case is at an angle of attack of 4° and therefore produces lift (C^ = 0.43). 

The airfoil is again almost completely immersed in supersonic flow, which 
begins at about 5% of chord (upper surface) and continues through oblique 
shocks at the trailing edge. A more complete description of the flow field is 
shown by the Mach number contour map displayed in figure 8. The oblique shock 
emanating from the trailing edge upper surface has been strengthened by the 
addition of circulation while the oblique shock emanating from the trailing 
edge lower surface has weakened and is almost nonexistent. The normal shock 
above the airfoil is much stronger than the normal shock below the airfoil. 
Except for these changes, the basic structure of the fishtail shock pattern is 
the same as for the first case. 

This is the first time in transonic flow computations that calculations 
such as those of the last two cases have been computed using an exact airfoil 
mapping. Without the newly developed rotated difference scheme these calcula- 
tions would have been unstable. The difficulty of these cases demonstrates 
the reliability of the present transonic flow solution procedure. 

Supersonic free-stream solutions— The pressure coefficient distribution 
for the last case (Moo = 1.15, a = 0°) is shown in figure 9. This case like 
the first is nonlifting and has a large amount of wave drag (Cp = 0.931). The 
decrease in wave drag for this case relative to the first fits in well with 
the established Cp vs M trend. An oscillation in the pressure distribution 
at the trailing edge is exhibited for this case and is attributed to the inter- 
action of the shock wave and the mapping singularity at the airfoil trailing 
edge. A more complete description of the entire flow field is shown by the 
Mach number contour map displayed in figure 10. A detached bow shock wave, 
characteristic of supersonic free streams, exists just upstream of the airfoil. 
Weak oblique shocks emanate from both the upper and lower surfaces at the air- 
foil trailing edge. The flow is supersonic at approximately 85% of all grid 
points in the flow field. The region of flow just downstream of the normal 
part of the bow shock and in the vicinity of the leading edge stagnation point 
is the only subsonic region. Thus, the finite-difference scheme over almost 
the entire flow field is first-order accurate. 

Some difficulties were experienced with this calculation and are attrib- 
uted to the use of standard subsonic boundary conditions on the outer boundary 
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(i.e., explicit specification of the fr 3. -stream solution). For instance, 
unstable results were obtained with the standard mesh. To stabilize this cal- 
culation the outer boundary had to be moved to a radius of 12 chords away from 
the airfoil instead of the usual 6 chords. In addition, the remesh option, 
which provides special grid point clustering at the airfoil (described in 
ref. 4), had to be used. Because supersonic free-stream cases are not of 
primary importance in the present study, no attempts have been made to stabi- 
lize this calculation on the standard mesh with a more appropriate supersonic 
outer boundary condition. 


C. Convergence Speed 

In the last section several classically difficult test cases have been 
presented in an attempt to establish the high level of reliability associated 
with the present new algorithm. Virtually any airfoil problem in the tran- 
sonic range (i.e., subsonic free stream) is deemed solvable by the present 
technique. In this section the convergence rate of these solutions is exam- 
ined. This topic was the subject of discussion in reference 4, where only 
moderately difficult cases were considered. An example of a particularly 
rapid convergence for an "easy" calculation is shown in figures 11 and 12. 

(The residual history for this calculation is given in ref. 4.) This calcula- 
tion is subcritical (Moo = 0.72) and nonlifting (a = 0°, NACA 0012). The pres- 
sure coefficient distributions after 9 and 15 iterations, with only every other 
point plotted for clarity, are compared with the solidly converged solution in 
figures 11 and 12, respectively. The solution trend is established after only 
9 iterations, and plottable accuracy is established after only 15 iterations. 
The latter solution corresponds to a reduction in the maximum residual of 
about one and a half orders of magnitude and represents only about 2 sec of 
computer time. The convergence speed exhibited in this calculation represents 
the ultimate speed achievable with the present algorithm. 

Of course, this type of convergence is not obtained in most cases. Both 
the existence of lift and supersonic regions of flow slow down convergence. 

An example of this is the convergence rate of the first case considered, i.e., 
an NACA 0012 airfoil, = 0.75, and a = 2° (see section III. A). Convergence 
histories for this case are displayed in figure 13. The lift coefficient (C^) , 
number of supersonic points (NSP) , and maximum residual (l^lmax^* normalized 
to 1.0 for the first iteration, are all plotted vs iteration number (n) . 

These curves, as well as all convergence history curves presented in this sec- 
tion, are constructed by plotting the appropriate quantity after every eighth 
iteration. The convergence rate in all cases has been approximately optimized 
by a trial and error process. Both NSP and Cj^ climb rapidly and reach 1% 
of their final values within 48 iterations. The lift coefficient overshoots 
slightly, which is characteristic of both Cj^ and NSP for many calculations. 
The flow field is essentially converged at this point even though the maximum 
residual has dropped by less than an order of magnitude. 

Establishment of crude convergence for such a small reduction in the max- 
imum residual is a characteristic behavior of many calculations using the 
present algorithm, especially for cases involving a large amount of supersonic 
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flow. For example, convergence history curves for the first two cases dis- 
cussed in section III.B (fishtail shock cases) are presented in figures 14 
and 15. The maximum residual in each case does not even start to i_i'op until 
after the solutions are approximately converged, l ..ring the initial phase of 
convergence in which the residual does not drop, the shock sonic line's posi- 
tion is rapidly being adjusted. This excites high-frequency errors, and there- 
fore, keeps the residual artificially high, even though the error is being 
reduced. The convergence history displayed in figure 14 indicates convergence 
in approximately 80 iterations, which is signified by the constant value of 
NSP. At this point the residual starts dropping exponentially, taking approx- 
imately 220 iterations to reach a six-order-of-magnitude reduction. The con- 
vergence history displayed in figure 15 indicates convergence in approximately 
100 iterations, and is signified by constant values of both NSP and 
Shortly after this point the residual starts dropping exponentially, ;aking 
approximately 300 iterations to reach a six-order-of-magnitude reduction in 
maximum residual. The degradation in convergence speed of the latter case is 
associated with the addition of lift into the problem, and more specifically, 
by the approximate technique for updating the circulation in the vortex solu- 
tion (e.g., see ref. 4). For larger amounts of lift the convergence rate will 
be even slower. 

The last convergence history curve presented in figure 16 is associated 
with the last test case discussed in section III.B (NACA 0012, = 1.15, 

a = 0°). Because for this case all grid points in the flow field are initially 
supersonic, the value of the NSP parameter decreases (instead of increasing) 
to the correct value as the solution converges. This solution is essentially 
converged in 100 iterations. The maximum residual again hesitates in the 
initial phase of convergence and then drops exponentially, taking approxi- 
mately 330 iterations to reach a six-order-of-magnitude reduction. 


IV. CONCLUSIONS 


A fast, implicit algorithm for solving the conservative full potential 
equation with effective rotated differencing has been presented. The rotated 
differencing is simply achieved by an upwind evaluation of the density coeffi- 
cient along both coordinate directions. This provides an effective upwind 
difference of the streamwise terms for any orientation of the velocity vector 
(i.e., rotated differencing), and, thereby, greatly enhances the reliability 
of the present algorithm. Use of the newly developed rotated differencing 
scheme has been instrumental in computing a number of classically difficult 
test cases including several cases with "fishtail" shock-wave patterns and a 
case with a detached bow shock wave. This represents the first time such cal- 
culations have been computed using the conservative full potential equation 
with an exact airfoil mapping. 

Results indicate that even for the classically difficult test cases pre- 
sented herein, the solution convergence rate is maintained at an exceptionally 
high level, which is indicative of the present fully implicit approximate fac- 
torization (AF2) algorithm. A characteristic of the AF2 convergence history 


15 



for many cases containing a large amount of supersonic flow (i.e., strong 
shock waves), is that the solution is essentially converged before the maximum 
residual even begins to drop. This indicates that the high-frequency error 
content of the solution does not decay until the solution is essentially con- 
verged. Because of this feature, the need for a different convergence cri- 
teria, other than a specified reduction in the maximum residual, is indicated. 
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Figure 1.- Numerically generated transformation, (x,y) ^ (^,n). 
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Figure 2.- Detail of a typical finite-difference i - i (NACA 0012 airfoil). 
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Figure 3.- Pressure coefficient distribution, artificial viscosity version 1, 

NACA 0012 airfoil, = 0.75, a = 2°. 
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Figure 12.- Intermediate solution comparison, NACA 0012 airfoil, = 0.72, 

a = 0°. 
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Figure 14.- Convergence history, NACA 0012 airfoil, = 0.98, a = 0°. 
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Figure 16.- Convergence history, I 
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